Add "3D Auto-correlograms" (3D ACGs) postprocessing extension by fededagos · Pull Request #3860 · SpikeInterface/spikeinterface

fededagos · 2025-04-11T09:09:09Z

Following conversations at Cosyne 2025 with @alejoe91, here is a first potential implementation of the "3D autocorrelograms" as an analyzer postprocessing extension.

First introduced in Beau et al., 2025, 3D ACGs provide a rich representation of a neuron's firing rate which accounts for different firing modes.
They were also recently used by Yu et al., (2025) again as a powerful representation of firing rates that can be used e.g. by representation learning methods, the same way that Beau et al. did.

The implementation is adapted directly from NeuroPyxels (mainly maintained by @m-beau), which contains the original Python version.

The original conceptualization is from @herzfeldd, who I am also tagging here for visibility.

alejoe91 · 2025-04-11T10:51:19Z

-
-
-from spikeinterface.core.waveforms_extractor_backwards_compatibility import MockWaveformExtractor
+from joblib import Parallel, delayed


@fededagos joblib is not a requirement. Instead, we normally use the built-in ProcessPoolExecutor

Is it ok if I push to this PR directly to refactor the parallelization?

Yes that is okay! Thanks

Done :) You can check the changes if you're interested!

…into feature/3d_acg

zm711

A couple initial comments. I guess one other approach would be to try to numba this instead for a speed gain since we already use numba for the standard correlograms. I saw that you pushed changes Alessio while I was typing this so if you changed something I marked feel free to ignore.

zm711 · 2025-04-11T15:37:08Z

+    n_jobs : int, default: -1.
+        Number of parallel jobs to spawn to compute the 3D-ACGS on different units.


Suggested change

n_jobs : int, default: -1.

Number of parallel jobs to spawn to compute the 3D-ACGS on different units.

If we want the job kwargs we need to do the {} trick and pull the shared kwarg docs. We can help with that.

zm711 · 2025-04-11T15:39:02Z

+
+    # pre-compute time bins
+    winsize_bins = 2 * int(0.5 * window_ms * 1.0 / bin_ms) + 1
+    bin_times_ms = np.linspace(-window_ms / 2, window_ms / 2, num=winsize_bins)


Suggested change

bin_times_ms = np.linspace(-window_ms / 2, window_ms / 2, num=winsize_bins)

bin_times_ms = np.linspace(-window_ms / 2, window_ms / 2, num=winsize_bins)

should we automatically make that the user can't specify a binsize that is inappropriate. So if I give a bin_size of like 0.0003 ms then most recording equipment don't have that resolution.

zm711 · 2025-04-11T16:54:06Z

+                acgs_3d[unit_index, :, :] = acg_3d
+                firing_quantiles[unit_index, :] = firing_quantile
+    else:
+        # Process units in serial


Suggested change

# Process units in serial

# Process units serially

:)

zm711 · 2025-04-11T16:55:14Z

+    assert winsize_bins >= 1
+    assert winsize_bins % 2 == 1


We should have some error output for this case for ease of debugging in the future.

zm711 · 2025-04-11T16:57:43Z

+try:
+    import npyx
+
+    HAVE_NPYX = True
+except ModuleNotFoundError:
+    HAVE_NPYX = False
+


importlib.util.find_spec

But I haven't fixed testing yet so I can go back and fix that later to try to gain a tiny speed boost in testing.

samuelgarcia · 2025-06-02T07:55:36Z

Hi.
Thanks a lot for this implementation.
The varaible naming and code style is really in the mood of spikeinterface and super easy to read.
Thank you.

Any idea why Maxime Beau chas choosen th name 3D-ACG ?
The name is not a super intuitive name regarding what the method does no ?
But keeping the original name is good idea if it is already adopted by the community.

I also have the intuition that a pure python loop will be a bit slow and numba would help a lot here.
We can do it later.

Maybe we could implement also in the same PR a very simple matplotlib widget no ? like plot_3dacg.

I will make some comments directly in the code.

samuelgarcia · 2025-06-02T07:57:10Z

+    smoothing_factor: int = 250,
+    **job_kwargs,
+):
+    """


If this is internal I would not put long doc here.
Or maybe I would rename this compute_acgs_3d_sorting and expose it to the API.

samuelgarcia · 2025-06-02T07:59:54Z

+            kernel_size = int(np.ceil(smoothing_factor / bin_size))
+            half_kernel_size = kernel_size // 2
+            kernel = np.ones(kernel_size) / kernel_size
+            smoothed_firing_rate = np.convolve(firing_rate, kernel, mode="same")


Maybe scipy.signal.oaconvolve could be faster here.

samuelgarcia · 2025-06-02T08:00:10Z

+        if isinstance(smoothing_factor, (int, float)) and smoothing_factor > 0:
+            kernel_size = int(np.ceil(smoothing_factor / bin_size))
+            half_kernel_size = kernel_size // 2
+            kernel = np.ones(kernel_size) / kernel_size


Why not precomputing the kernel before the unit_id loop ?

samuelgarcia · 2025-06-02T08:00:27Z

+        if isinstance(smoothing_factor, (int, float)) and smoothing_factor > 0:
+            kernel_size = int(np.ceil(smoothing_factor / bin_size))
+            half_kernel_size = kernel_size // 2
+            kernel = np.ones(kernel_size) / kernel_size


Maybe we could also have a gaussian kernel no ?

samuelgarcia · 2025-06-02T08:02:13Z

+            # Correct manually for possible artefacts at the edges
+            for i in range(kernel_size):
+                start = max(0, i - half_kernel_size)
+                stop = min(len(firing_rate), i + half_kernel_size)
+                smoothed_firing_rate[i] = np.mean(firing_rate[start:stop])
+            for i in range(len(firing_rate) - kernel_size, len(firing_rate)):
+                start = max(0, i - half_kernel_size)
+                stop = min(len(firing_rate), i + half_kernel_size)
+
+                smoothed_firing_rate[i] = np.mean(firing_rate[start:stop])


OK with this but this prevent using other kernel for smoothing.
Is this border handling usefull ?

m-beau · 2025-06-08T21:23:53Z

Hey there, thanks for the input Samuel!

David Herzfeld called it '3D-ACG' because regular ACGs are 2-dimensional (time-shift, correlation/rate), and those are simply stratified by instantaneous firing rate. Thus they're graphically plotted as heatmaps, which intuitively relate to 3D arrays!

Because that's what we called them in Beau et al. 2025, and also in David's papers, I'm reluctant to rename them :)

Super happy to see this merged into spike-interface! By the way, @alejoe91 we should get together with Julie Fabre and chat about quality metrics / bombcell but also about a potential Phy3 someday soon. Cyrille Rossant told me that you were interested in writing it up - we had something similar in mind, we should merge forces!

Maxime

herzfeldd · 2025-06-08T21:38:08Z

Thanks, Maxime!

I'm reluctant to rename them as well, as they are now '3D-ACGs' in multiple papers.

I will note while they are a creation of myself and Nate Hall, they were officially named by Steve Lisberger. I initially called them firing rate versus correlogram plots (because they trivially extend to CCGs as well as ACGs) but 3D-ACG is a lot fewer words...

alejoe91 · 2025-06-09T11:18:01Z

Super happy to see this merged into spike-interface! By the way, @alejoe91 we should get together with Julie Fabre and chat about quality metrics / bombcell but also about a potential Phy3 someday soon. Cyrille Rossant told me that you were interested in writing it up - we had something similar in mind, we should merge forces!

Of course! I'll ping you on slack :)

fededagos added 2 commits April 11, 2025 10:45

Added 3D ACGs, references, tests

d78573c

Merging changes to main branch

b29c1c5

alejoe91 reviewed Apr 11, 2025

View reviewed changes

alejoe91 added the postprocessing Related to postprocessing module label Apr 11, 2025

alejoe91 added 2 commits April 11, 2025 17:29

Use SI parallel machinery

5bb0460

Merge branch 'main' into feature/3d_acg

c1b47f1

alejoe91 requested a review from zm711 April 11, 2025 15:30

alejoe91 added 2 commits April 11, 2025 18:18

Fix tests

d3fddee

Merge branch 'feature/3d_acg' of github.com:fededagos/spikeinterface …

81b7373

…into feature/3d_acg

zm711 reviewed Apr 11, 2025

View reviewed changes

Assertion messages, importlib, and job_kwargs docs

55adef6

samuelgarcia reviewed Jun 2, 2025

View reviewed changes

alejoe91 added this to the 0.103.0 milestone Jun 11, 2025

fix conflicts

18bb2f7

samuelgarcia approved these changes Jun 12, 2025

View reviewed changes

fix conflicts 2

b15411c

alejoe91 merged commit 8dfef70 into SpikeInterface:main Jun 12, 2025
15 checks passed



		from spikeinterface.core.waveforms_extractor_backwards_compatibility import MockWaveformExtractor
		from joblib import Parallel, delayed

		n_jobs : int, default: -1.
		Number of parallel jobs to spawn to compute the 3D-ACGS on different units.

	bin_times_ms = np.linspace(-window_ms / 2, window_ms / 2, num=winsize_bins)
	bin_times_ms = np.linspace(-window_ms / 2, window_ms / 2, num=winsize_bins)

Conversation

fededagos commented Apr 11, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

zm711 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

samuelgarcia commented Jun 2, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

m-beau commented Jun 8, 2025

Uh oh!

herzfeldd commented Jun 8, 2025

Uh oh!

alejoe91 commented Jun 9, 2025

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

6 participants